Syntax

A quick recap of the syntax used within formulas. The ~ (tilde) sign generally differentiates between dependent and independent variables. The + sign is used to distinguish between independent variables. Here are a few more common formula inputs that are useful to know:

y ~ x1 + x2 # standard formula

y ~ . # include all other columns in data as x vars

y ~ x1 * x2 # interaction terms

y ~ x + I(x^2) # higher order terms

Formulas are generally straightforward to work with. To avoid any potential errors, you should follow these two steps:

Step 1: Create a string containing the written formula.

formula_string <- paste("life_expectancy", "~ gdp")
formula_string

## [1] "life_expectancy ~ gdp"

Step 2: Transform the string into R’s formula class. Only then will your modeling function accept it as input.

form <- as.formula(formula_string) # needs to be transformed into formula class
form

## life_expectancy ~ gdp

Here is the proof that both ways of doing it are identical.

reg1 <- lm(form, data = life_expec_dat) 
reg2 <- lm(life_expectancy ~ gdp, data = life_expec_dat)

reg1$coefficients == reg2$coefficients

## (Intercept)         gdp 
##        TRUE        TRUE

Iteration

In instances where you might need to fit many models, formulas can help with making this process iterative.

Often, the modelling process requires you to run the same specification with multiple configurations of both dependent and independent variables. Model formulas make running many similar models much easier:

Step 1: Define a function that let’s you plug in different variables for x or y. The function should take the x and/or y variable(s) as a string and return a model object.

# function to include different y variables
lm_fun_iter_y <-  function(dep_var) {
  lm(as.formula(paste(dep_var, "~ gdp")), data = life_expec_dat)
}

# function to include different x variables 
lm_fun_iter_x <- function(indep_var) {
  lm(as.formula(paste("life_expectancy ~", paste(indep_var, collapse = "+"))), data = life_expec_dat)
}

Notice: Considering your model will likely include many independent variables, it’s unlikely that you’ll run a simple bivariate regression. Therefore, we need to use a nested paste function that combines (with collapse()) all independent variables contained in a character vector with +.

Step 2: Use purrr::map() (refer to lab session 4) to iteratively apply the model to each variable contained in a character vector of variables:

# create vector of variables to iterate over
vars <- life_expec_dat |> 
  dplyr::select(-c("life_expectancy", "country")) |> 
  names()

vars

##  [1] "year"                            "status"                         
##  [3] "adult_mortality"                 "infant_deaths"                  
##  [5] "alcohol"                         "percentage_expenditure"         
##  [7] "hepatitis_b"                     "measles"                        
##  [9] "bmi"                             "under_five_deaths"              
## [11] "polio"                           "total_expenditure"              
## [13] "diphtheria"                      "hiv_aids"                       
## [15] "gdp"                             "population"                     
## [17] "thinness_1_19_years"             "thinness_5_9_years"             
## [19] "income_composition_of_resources" "schooling"                      
## [21] "large"

# run a bivariate model for each column
biv_model_out <-  vars |>
  purrr::map(lm_fun_iter_x)

biv_model_out

## [[1]]
## 
## Call:
## lm(formula = as.formula(paste("life_expectancy ~", paste(indep_var, 
##     collapse = "+"))), data = life_expec_dat)
## 
## Coefficients:
## (Intercept)         year  
##    -635.872        0.351  
## 
## 
## [[2]]
## 
## Call:
## lm(formula = as.formula(paste("life_expectancy ~", paste(indep_var, 
##     collapse = "+"))), data = life_expec_dat)
## 
## Coefficients:
##      (Intercept)  statusDeveloping  
##             79.2             -12.1  
....

This returns output from all of the models that you’ve just run. Given that there are now quite a few models, it might be a good idea to keep the names of the independent variables that you feed into the model. You can do this with purrr::set_names().

# run a bivariate model for each column
biv_model_out_w_names <-  vars |>
  purrr::set_names() |> 
  purrr::map(lm_fun_iter_x)

biv_model_out_w_names$year # we can now easily index independent vars

## 
## Call:
## lm(formula = as.formula(paste("life_expectancy ~", paste(indep_var, 
##     collapse = "+"))), data = life_expec_dat)
## 
## Coefficients:
## (Intercept)         year  
##    -635.872        0.351

Say that you’re interested in seeing what effect gdp has on life_expectancy, but you suspect that other covariates such as alcohol consumption, health care expenditure, average body mass index and the propagation of aids also have an effect.

On top of this, it’s possible that these variables might have different effects on life_expectancy depending on how they interact with each other. We can determine which independent variables are robustly associated with the dependent variable by generating models for every possible combination of independent variables and then looking into the distribution of effects across these models.

Let’s review the code introduced in the lecture:

Step 1: Find out how many possible combinations of independent and dependent variables exist:

combinations <- 
  purrr::map(1:5, function(x) {combn(1:5, x)}) |> # create all possible combinations (assume we have 5 covariates)
  purrr::map(ncol) |> # extract number of combinations
  unlist() |> # pull out of list structure
  sum() # compute sum

combinations

## [1] 31

Step 2: Write a function that can run all possible combinations in one go:

combn_models <- function(depvar, covars, data) {
  
  # initialize empty list
  combn_list <- list()
  
  # generate list of covariate combinations
  for (i in seq_along(covars)) {
    combn_list[[i]] <- combn(covars, i, simplify = FALSE)
  } # what is combn?
  
  # unlist to make object easier to work with
  combn_list <-
    unlist(combn_list, recursive = FALSE) 
  
  # function to generate formulas
  gen_formula <- function(covars, depvar) {
    form <-
      as.formula(paste0(depvar, " ~ ", paste0(covars, collapse = "+")))
    form
  }
  
  # generate formulas
  formulas_list <-
    purrr::map(combn_list, gen_formula, depvar = depvar)
  
  # run models
  models_list <- purrr::map(formulas_list, lm, data = data)
  models_list
}

Step 3: Apply the function to our case:

depvar <- "life_expectancy"
covars <-  c("gdp", "percentage_expenditure", "alcohol", "bmi", "hiv_aids")
multiv_model_out <- combn_models(depvar = depvar, covars = covars, data = life_expec_dat)

# check whether we have the correct number of models
length(multiv_model_out)

## [1] 31

As you can see, the function introduced in the lecture is quite powerful, but it is by no means the only option for testing multiple model specifications. Depending on the type of analysis that you plan to perform, certain R packages enable you run many different combinations of regression models with just a few lines of code. For example, the we can perform multiple estimations with the fixest package.

Tidying Model Objects

Luckily, we don’t have to deal with this. We can use the tidyverse’s broom package which is part of a collection of packages used for modeling that adhere to tidyverse principles referred to as tidymodels. Tidymodels offers a lot of cool packages that have their own tutorials, so check them out! 💪.

In today’s session, we will only focus on broom. As a refresher, here are the three key broom functions that you should know:

1. broom::tidy(): Summarizes information about model components.

broom::tidy(multiv_model_out[[3]], conf.int = TRUE, conf.level = 0.95)

## # A tibble: 2 × 7
##   term        estimate std.error statistic   p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)   64.8      0.253      256.  0           64.3       65.3 
## 2 alcohol        0.955    0.0412      23.1 2.11e-108    0.874      1.04

2. broom::glance(): Reports information about the entire model.

broom::glance(multiv_model_out[[3]])

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic   p.value    df logLik    AIC    BIC
##       <dbl>         <dbl> <dbl>     <dbl>     <dbl> <dbl>  <dbl>  <dbl>  <dbl>
## 1     0.164         0.164  8.73      536. 2.11e-108     1 -9807. 19621. 19639.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

3. augment(): Adds information about observations to a dataset.

broom::augment(multiv_model_out[[3]], se_fit = TRUE) 

## # A tibble: 2,735 × 10
##    .rownames life_expectancy alcohol .fitted .se.fit .resid     .hat .sigma
##    <chr>               <dbl>   <dbl>   <dbl>   <dbl>  <dbl>    <dbl>  <dbl>
##  1 1                    65      0.01    64.8   0.253  0.227 0.000838   8.74
##  2 2                    59.9    0.01    64.8   0.253 -4.87  0.000838   8.74
##  3 3                    59.9    0.01    64.8   0.253 -4.87  0.000838   8.74
##  4 4                    59.5    0.01    64.8   0.253 -5.27  0.000838   8.74
##  5 5                    59.2    0.01    64.8   0.253 -5.57  0.000838   8.74
##  6 6                    58.8    0.01    64.8   0.253 -5.97  0.000838   8.74
##  7 7                    58.6    0.01    64.8   0.253 -6.17  0.000838   8.74
##  8 8                    58.1    0.03    64.8   0.252 -6.69  0.000834   8.74
##  9 9                    57.5    0.02    64.8   0.253 -7.28  0.000836   8.74
## 10 10                   57.3    0.03    64.8   0.252 -7.49  0.000834   8.74
## # ℹ 2,725 more rows
## # ℹ 2 more variables: .cooksd <dbl>, .std.resid <dbl>

Doesn’t this look so much tidier?

Broom with Many Models

Where broom really shines is when it comes to dealing with multiple models like we have in our multiv_model_out. It allows us to easily move away from hard-to-deal-with lists while retaining each of the three important model outputs. Since our example replicates several variables in different models, we include the .id argument in purrr::map(). This adds a column to our tibble which will contain a model identifier. In our example, the identifier of each model corresponds with the index number of the model output.

multiv_model_out_broom <- purrr::map_dfr(multiv_model_out[1:3], # the model objects
                                         broom::tidy, # the function to iterate with
                                         .id = "model_type") # id column

multiv_model_out_broom

## # A tibble: 6 × 6
##   model_type term                    estimate std.error statistic   p.value
##   <chr>      <chr>                      <dbl>     <dbl>     <dbl>     <dbl>
## 1 1          (Intercept)            67.0      0.194         346.  0        
## 2 1          gdp                     0.000312 0.0000120      25.9 2.71e-131
## 3 2          (Intercept)            67.9      0.174         391.  0        
## 4 2          percentage_expenditure  0.00183  0.0000817      22.3 2.77e-102
## 5 3          (Intercept)            64.8      0.253         256.  0        
## 6 3          alcohol                 0.955    0.0412         23.1 2.11e-108

Now that we have multiple model outputs in a format that we can work with, it is time to start thinking about visualizing our results. You can plot the residuals, the goodness of fit statistics etc. Similarly, you could also use broom::glance() to provide some robustness statistics or descriptive tables.

The next section will focus on how to present your results!

Using Regression Tables

The workhorse function in modelsummary is unsurprisingly modelsummary(). 🤔

Here we only need to input the models list as the function does all of the tidying for us.

modelsummary::modelsummary(multiv_model_out[[31]], 
                           output = "kableExtra") # can take multiple formats (e.g., .tex, .doc, .pdf)

It looks pretty solid already, but can definitely be improved upon.

For example, we can set acceptable number of digits:

modelsummary::modelsummary(multiv_model_out[[31]], 
                           output = "kableExtra",
                           fmt = "%.2f") # 2-digits and trailing 0

We can also report only the C.I. by getting rid of the p-values and intercept term:

model_table <- modelsummary::modelsummary(multiv_model_out[[31]], 
                                          output = "kableExtra",
                                          fmt = "%.2f",  # 2-digits and trailing 0  
                                          estimate  = "{estimate}", 
                                          statistic = "conf.int",
                                          coef_omit = "Intercept") 

Wouldn’t it look more professional if our coefficients began with capital letter? And how about we get rid of some of the bloated goodness of fit statistics?

mod_table <- modelsummary::modelsummary(multiv_model_out[31], 
                                        output = "default",
                                        fmt = "%.2f",  # 2-digits and trailing 0  
                                        estimate  = "{estimate}",
                                        statistic = "conf.int",
                                        coef_omit = "Intercept",
                                        coef_rename=c("gdp"="Gdp", 
                                                      "bmi"="Avg. BMI", 
                                                      "alcohol" = "Alcohol Consum.",
                                                      "hiv_aids"= "HIV cases", 
                                                      "percentage_expenditure" = "Health Expenditure (% of GDP)",
                                                      "life_expectancy" = "Life Expectancy"),
                                        gof_omit = 'DF|Deviance|Log.Lik|AIC|BIC',
                                        title = 'An Improved Regression Table')

mod_table

When using the kableExtra package, you can even post-process your table:

mod_table |> 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F, fixed_thead = T) |> 
  kableExtra::row_spec(3, color = 'red') |>
  kableExtra::row_spec(5, background = 'lightblue')

Using Plots 🎨

If you are looking for an alternative way to present your results, you might also consider a nice coefficient plot.

Again, modelsummary provides a very accessible function for this: modelplot(). As with the regression table, you input the Model List and not the tidied output. The tidying happens under the hood.

# to avoid overplotting, we define a smaller list of models
models <- list(
  "Bivariate OLS" = lm(life_expectancy ~ gdp, data = life_expec_dat),
  "Multivariate OLS" = lm(life_expectancy ~ gdp + total_expenditure + population + alcohol + bmi, data = life_expec_dat),
  "Interaction" = lm(life_expectancy ~  gdp + total_expenditure + population + gdp*total_expenditure, data = life_expec_dat),
  "Polynomial" = lm(life_expectancy ~ gdp + I(gdp^2), data = life_expec_dat)
)

# plot 
modelsummary::modelplot(models, coef_omit = "Intercept") 

Now we just need to customize the plot to our liking:

modelsummary::modelplot(models, coef_omit = "Intercept") +
  labs(x = 'Coefficients', 
       y = 'Term names',
       title = 'Linear regression models of "Life expectancy"',
       caption = "Comparing multiple models, the data originated at the WHO") 

We will cover data visualization with R in much greater detail in the next session, so hopefully the section on coefficient plot will be made clearer then!

Crafting formulas

The basic structure of a formula:

y ~ x

Running a model with a formula is straightforward. Note that we don’t even have to put the formula in parentheses:

lm(arr_delay ~ distance + origin, data = flights)

## 
## Call:
## lm(formula = arr_delay ~ distance + origin, data = flights)
## 
## Coefficients:
## (Intercept)     distance    originJFK    originLGA  
##    13.41405     -0.00405     -2.70426     -4.45617

A more explicit way to write formulas:

fmla <- as.formula("arr_delay ~ distance + origin")
class(fmla)

## [1] "formula"

lm(fmla, data = flights)

## 
## Call:
## lm(formula = fmla, data = flights)
## 
## Coefficients:
## (Intercept)     distance    originJFK    originLGA  
##    13.41405     -0.00405     -2.70426     -4.45617

More formula syntax:

y ~ x1 + x2

## y ~ x1 + x2

y ~ x1 - x2 # ignore x2 in analysis

## y ~ x1 - x2

y ~ . # select all other variables in model matrix

## y ~ .

Interactions:

y ~ x1 * x2 # That's equivalent to y ~ x1 + x2 + x1*x2

## y ~ x1 * x2

y ~ x1:x2  # only interaction, no main effects

## y ~ x1:x2

As-is variables:

y ~ x + I(x^2)

## y ~ x + I(x^2)

Extract all variable names from a formula object:

all.vars(fmla)

## [1] "arr_delay" "distance"  "origin"

Incrementally modify/update a formula:

#?update

Advanced crafting of formulas: fitting multiple models

Create a formula for a model with a large number of variables:

xnam <- paste0("x", 1:25)
(fmla <- as.formula(paste("y ~ ", paste(xnam, collapse= "+"))))

## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + 
##     x12 + x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + 
##     x22 + x23 + x24 + x25

depvar <- "arr_delay"
covars <- c("dep_delay", "carrier", "origin", "air_time", "distance", "hour")

combn_models <- function(depvar, covars, data)
{
  combn_list <- list()
  # generate list of covariate combinations
  for (i in seq_along(covars)) {
    combn_list[[i]] <- combn(covars, i, simplify = FALSE)
  }
  combn_list <- unlist(combn_list, recursive = FALSE)
  # function to generate formulas
  gen_formula <- function(covars, depvar) {
    form <- as.formula(paste0(depvar, " ~ ", paste0(covars, collapse = "+")))
    form
  }
  # generate formulas
  formulas_list <- purrr::map(combn_list, gen_formula, depvar = depvar)
  # run models
  models_list <- purrr::map(formulas_list, lm, data = data)
  models_list
}

models_list <- combn_models(depvar = depvar,  covars = covars, data = flights)
length(models_list)

## [1] 63

Working with tidy model outputs

Inspecting models in R is straightforward with the summary() function.

lm(arr_delay ~ distance + origin, data = flights) |> summary()

## 
## Call:
## lm(formula = arr_delay ~ distance + origin, data = flights)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -89.0  -24.0  -11.8    7.3 1281.4 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.41405    0.17481    76.7   <2e-16 ***
## distance    -0.00404    0.00011   -36.9   <2e-16 ***
## originJFK   -2.70425    0.18871   -14.3   <2e-16 ***
## originLGA   -4.45617    0.19351   -23.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 44 on 327342 degrees of freedom
##   (9430 observations deleted due to missingness)
## Multiple R-squared:  0.0055, Adjusted R-squared:  0.00549 
## F-statistic:  604 on 3 and 327342 DF,  p-value: <2e-16

But often you want to post-process estimation results. So let’s examine the output of a model function more closely.

model_out <- lm(arr_delay ~ distance + origin, data = flights)
#class(model_out)
#str(model_out)

Ugh, that’s a lot of information in a list structure. Let’s use some helper functions to unpack these:

#coef(model_out)
#fitted.values(model_out)
#residuals(model_out)
#model.matrix(model_out)

We can also dig in the summary of the model.

#str(summary(model_out))
summary(model_out)$coefficients

##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   13.414    0.17481      77  0.0e+00
## distance      -0.004    0.00011     -37 5.5e-297
## originJFK     -2.704    0.18871     -14  1.5e-46
## originLGA     -4.456    0.19351     -23 3.0e-117

model_out <- lm(arr_delay ~ distance + origin, data = flights)
model_out

## 
## Call:
## lm(formula = arr_delay ~ distance + origin, data = flights)
## 
## Coefficients:
## (Intercept)     distance    originJFK    originLGA  
##    13.41405     -0.00405     -2.70426     -4.45617

summary(model_out)

## 
## Call:
## lm(formula = arr_delay ~ distance + origin, data = flights)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -89.0  -24.0  -11.8    7.3 1281.4 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.41405    0.17481    76.7   <2e-16 ***
## distance    -0.00404    0.00011   -36.9   <2e-16 ***
## originJFK   -2.70425    0.18871   -14.3   <2e-16 ***
## originLGA   -4.45617    0.19351   -23.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 44 on 327342 degrees of freedom
##   (9430 observations deleted due to missingness)
## Multiple R-squared:  0.0055, Adjusted R-squared:  0.00549 
## F-statistic:  604 on 3 and 327342 DF,  p-value: <2e-16

#?tidy.lm
broom::tidy(model_out)

## # A tibble: 4 × 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept) 13.4      0.175         76.7 0        
## 2 distance    -0.00405  0.000110     -36.9 5.53e-297
## 3 originJFK   -2.70     0.189        -14.3 1.46e- 46
## 4 originLGA   -4.46     0.194        -23.0 3.04e-117

broom::tidy(model_out, conf.int = TRUE, conf.level = 0.95)

## # A tibble: 4 × 7
##   term        estimate std.error statistic   p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept) 13.4      0.175         76.7 0         13.1      13.8    
## 2 distance    -0.00405  0.000110     -36.9 5.53e-297 -0.00426  -0.00383
## 3 originJFK   -2.70     0.189        -14.3 1.46e- 46 -3.07     -2.33   
## 4 originLGA   -4.46     0.194        -23.0 3.04e-117 -4.84     -4.08

#?glance.lm
broom::glance(model_out)

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic p.value    df    logLik     AIC    BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>     <dbl>   <dbl>  <dbl>
## 1   0.00550       0.00549  44.5      604.       0     3 -1706997.  3.41e6 3.41e6
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

#?augment.lm
broom::augment(model_out) |> slice(1:5)

## # A tibble: 5 × 10
##   .rownames arr_delay distance origin .fitted .resid       .hat .sigma   .cooksd
##   <chr>         <dbl>    <dbl> <chr>    <dbl>  <dbl>      <dbl>  <dbl>     <dbl>
## 1 1                11     1400 EWR       7.75   3.25 0.00000922   44.5   1.23e-8
## 2 2                20     1416 LGA       3.23  16.8  0.0000123    44.5   4.37e-7
## 3 3                33     1089 JFK       6.30  26.7  0.00000938   44.5   8.43e-7
## 4 4               -18     1576 JFK       4.33 -22.3  0.00000972   44.5   6.12e-7
## 5 5               -25      762 LGA       5.88 -30.9  0.00000989   44.5   1.19e-6
## # ℹ 1 more variable: .std.resid <dbl>

tt <- t.test(wt ~ am, mtcars)
tidy(tt)

## # A tibble: 1 × 10
##   estimate estimate1 estimate2 statistic    p.value parameter conf.low conf.high
##      <dbl>     <dbl>     <dbl>     <dbl>      <dbl>     <dbl>    <dbl>     <dbl>
## 1     1.36      3.77      2.41      5.49 0.00000627      29.2    0.853      1.86
## # ℹ 2 more variables: method <chr>, alternative <chr>

models_broom <- purrr::map(models_list, broom::tidy)
models_broom[[1]] # inspect one list entry

## # A tibble: 2 × 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)    -5.90  0.0330       -179.       0
## 2 dep_delay       1.02  0.000786     1296.       0

# turn into df by rbinding all tidied results
models_broom_df <- purrr::map_dfr(models_broom, rbind)

models_broom_df |> 
  dplyr::filter(!str_detect(term, "Intercept|carrier")) |>
  ggplot(aes(estimate)) + 
  geom_histogram(binwidth = .1) + 
  geom_vline(xintercept = 0, linetype="dashed") + 
  facet_grid(cols = vars(term), scales = "free_y") + 
  theme_minimal()

Reporting model results

model1_out <- lm(arr_delay ~ dep_delay, data = flights)
model2_out <- lm(arr_delay ~ distance, data = flights)
model3_out <- lm(arr_delay ~ dep_delay + distance, data = flights)

models <- list(model1_out, model2_out, model3_out)

modelsummary::modelsummary(models, 
                           estimate = "{estimate} [{conf.low}, {conf.high}]",
                           statistic = NULL,
                           gof_omit = ".+",
                           title = "Linear regression of flight delay at arrival (in mins)") 

cm <- c("distance" = "Distance",
        "originLGA" = "Origin: LGA",
        "originJFK" = "Origin: JFK")

modelsummary::modelplot(model_out, 
                        coef_omit = "Interc", 
                        coef_map = cm)

flights$distance_std <- effectsize::standardize(flights$distance)
model_out_std2 <-lm(arr_delay ~ distance_std + origin, data = flights)

cm <- c("distance_std" = "Distance",
        "originLGA" = "Origin: LGA",
        "originJFK" = "Origin: JFK")

modelsummary::modelplot(model_out_std2, coef_omit = "Interc", coef_map = cm) + 
  xlim(-5, .25) + 
  geom_vline(xintercept = 0, linetype="dashed") + 
  labs(title = "Linear regression of flight delay at arrival (in mins)",
       caption = "Data source: nycflights13 package") + 
  theme_minimal()

flights$distance1kmiles <- flights$distance / 1000
model_out_kmiles <-lm(arr_delay ~ distance1kmiles + origin, data = flights)
cm <- c("distance1kmiles" = "Distance (1k miles)",
        "originLGA" = "Origin: LGA",
        "originJFK" = "Origin: JFK")

modelsummary::modelplot(model_out_kmiles, coef_omit = "Interc", coef_map = cm) + 
  xlim(-5, .25) + 
  geom_vline(xintercept = 0, linetype="dashed") + 
  labs(title = "Linear regression of flight delay at arrival (in mins)",
       caption = "Data source: nycflights13 package") + 
  theme_minimal()

modelsummary::modelplot(models, coef_omit = "Interc") + 
  facet_grid(~model)

library(effectsize)
#?standardize
#model_out_std <- standardize_parameters(model_out)
#modelplot(model_out_std, coef_omit = "Interc", coef_map = cm) + 
#  labs(title = "Linear regression of flight delay at arrival (in mins)",
#       caption = "Data source: nycflights13 package") + 
#  theme_minimal()

model_out |> 
  broom::augment() |>
  dplyr::slice_sample(n = 10000) |> 
  ggplot(aes(x = .fitted, y = arr_delay)) +
  geom_point(alpha =  0.25) +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Fitted vs. true values, lm(arr_delay ~ distance + origin)",
       caption = "Data source: nycflights13 package") + 
  xlab("Fitted values") + 
  ylab("Arrival delay (in mins)") + 
  theme_minimal()

model_out_improve <- lm(arr_delay ~ dep_delay + distance + origin, data = flights)

model_out_improve |> 
  broom::augment() |>
  dplyr::slice_sample(n = 10000) |> 
  ggplot(aes(x = .fitted, y = arr_delay)) +
  geom_point(alpha =  0.25) +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Fitted vs. true values, lm(arr_delay ~ dep_delay + distance + origin)",
       caption = "Data source: nycflights13 package") + 
  xlab("Fitted values") + 
  ylab("Arrival delay (in mins)") + 
  theme_minimal()